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The Nose Hamiltonian is adapted, leading to a deriva- 
tion of the Nose-Hoover equations of motion which does not 
involve time transformations, and in which the degree of free- 
dom corresponding to the external reservoir is treated on the 
same footing as those of the rest of the system. In this form 
it is possible to prove the conjugate pairing rule for Lyapunov 
exponents of this system. 
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I. INTRODUCTION 



(5) 



The Nose Hamiltonian 



H N (q, s; 7r,p s ; A) = ^ ^3 + tp(q) + + gkTlns , 

(i) 

is used to model a system of N particles interacting with 
a thermal reservoir at temperature T, represented by 
the coordinate s with its conjugate momentum p s . g 
is a constant which depends on the number of degrees 
of freedom of the system. If we interpret the time vari- 
able A as a nonphysical parameter, with physical time 
intervals defined by dt = dX/s and physical momentum 
Pi = dqi/dt = Tti/s, a uniform (microcanonical) prob- 
ability measure on the full (q, 7r, s,p s ) phase space re- 
duces to a canonical probability measure on the system 
variables (q, p) , as long as g = 3N + 1 . Thus a molecu- 
lar dynamics simulation may be used to model canonical, 
as well as microcanonical ensembles, as long as the ex- 
tended system dynamics is ergodic. Refer to Ref. for a 
more complete discussion and review. Numerical results 
suggest that the ergodic assumption may be reasonable 
for all but the very smallest systems 0. It is possible 
make a number of modifications of the time and momen- 
tum variables in Hn to obtain other dynamical systems 
which generate the canonical ensemble |5]](| . 

Hoover pointed out that the equations in terms of 
the physical variables (q, p, t) and C = Ps/Q take on a 
particularly simple form: 



q, = Pi /mi 

p t = -Vnp - (Pi 

c = 




gkT 



(2) 
(3) 

(4) 



Note that the first three equations form a closed set; s is 
now redundant. In this form, it is apparent that £ acts 
as a kind of thermostat acting on the kinetic energy of 
the iV-particle system. £ is proportional to the difference 
between the kinetic energy and gkT/2, so that when the 
kinetic energy rises above this value, £ increases, raising 
the damping term in the equation for p, and thus reduc- 
ing the kinetic energy. In the Hoover representation the 
equilibrium distribution is changed because of the time 
scaling, however the canonical distribution in the system 
variables is recovered by setting g = 3N. 

The real utility of this approach, however, is to 
nonequilibrium systems. In a nonequilibrium simula- 
tion in which energy is being pumped through the sys- 
tem such a thermostat permits the system to approach a 
steady state while remaining homogeneous; there are al- 
ternative approaches involving boundaries which do not 
share this property [|s| — |l~0|] , Another important ther- 
mostatting method based on Gauss' principle of least 
constraint jll) uses £ as an explicit function of the coor- 
dinates rather than a variable in its own right j p2| . The 
Nose-Hoover and Gaussian thermostats give the same 
averages and time correlation functions in the thermo- 
dynamic limit |l3|| . The important point to note here is 
that any representation of a nonequilibrium steady state 
must contain some reference to an external heat reser- 
voir. A possible advantage of the Nose-Hoover scheme 
over the Gaussian approach is that this reservoir is ex- 
plictly included as a separate degree of freedom, so that it 
may be treated on a similar footing as the rest of the sys- 
tem. Recently a number of modifications and extensions 
to the Nose-Hoover method have been proposed [p^Tl7[ . 

Lyapunov exponents (defined in Sect. Ill) are impor- 
tant in the study of nonequilibrium systems, giving in- 
formation on the chaotic instability, and providing an 
important link between the microscopic and macroscopic 
properties, since the sum of the exponents gives the av- 
erage rate of phase space expansion, which can then be 
related to entropy production, and hence transport co- 
efficients. The conjugate pairing rule for Lyapunov ex- 
ponents is the property that there is a constant C such 
that for every exponent A, C— A is also an exponent, with 
the possible exception of one or two exponents which are 
fixed to be zero by symmetry considerations. It was first 
discussed in Refs. |l~2| , |l8| |. The conjugate pairing rule 
permits the sum of the exponents to be calculated from 
the largest and smallest, which are the easiest to eval- 
uate numerically. Hamiltonian systems obey conjugate 
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pairing with C — 0, that is, the exponents come in ± 
pairs |19f |. Systems with a constant damping factor pair 
with C proportional to this factor ^0|. Recently the 
conjugate pairing rule has been shown to hold for sys- 
tems containg a Gaussian thermostat M], where C is 
minus the average value of the thermostatting multiplier 
a (analogous to C here). In this case there are two zero 
exponents which do not pair, which arise from the time 
translation symmetry, and the conserved kinetic energy. 

Note that the Lyapunov exponents depend on the vari- 
ables used to define the phase space if the equations re- 
lating different coordinate systems involve exponential 
functions of time. This means that, although it is triv- 
ial to prove that the Lyapunov exponents obtained us- 
ing the original Nose variables pair to zero, because the 
equations of motion are derived from a Hamiltonian, it 
is much harder to make statements about the Hoover 
variables, particularly since a different time variable is 
used. Here the total phase space contraction, which is 
given by minus the sum of the Lyapunov exponents, is 
proportional to the average of which is nonzero for a 
nonequilibrium steady state. Since the sum of the expo- 
nents is less than zero, it is clear that the exponents are 
quite different to the Nose values, for which the sum is 
trivially zero. 

Sect. H of this paper shows how to write the Nose 
Hamiltonian in a form which treats the system and reser- 
voir variables on equal footing, leading to a proof of the 
conjugate pai ring rule in Sect. HI , along the lines of the 
proof in Ref. || . 



II. UNIFIED HAMILTONIAN FORMALISM 

In the form given by Nose (Eq. [TJ) , the reservoir variable 
s is treated quite differently to the other coordinates in 
that the system kinetic term is divided by s 2 wheras the 
reservoir kinetic term is not. In addition, the time A 
corresponding to the Hamiltonian does not correspond 
to physical time. In this section Hn is transformed to 
alleviate both of these deficiencies. 

The unified form of the Hamiltonian is obtained by 
transforming to a new coordinate a = Ins. This type of 
transformation is described in Sect. 9.2 of Ref. plfl , and 
uses a generating function of the form 



N 



F 2 (c[,s;Tv,p a ;X) = • TVi+p a lns 



(6) 



leading to a new momentum p a — sp s , and transformed 
Hamiltonian 



H T {q,(7]TT,p a ;X) 



2 [fern, Q, 
<p(q) + gkTa. 



(7) 



Note that the form of the potential is also simpler in this 
representation. The masses may also be scaled out: Con- 
struct 3N + 1 dimensional vectors X = {ciiy/ml, &\/~Q) 
and P = (iTi/ \fmliPa / \/Q), and write ^(X) = — a and 
0(X) = ip(q) + gkTa. The Hamiltonian may then be 
written in a unified form as 



i^(X;P;A) = ^e 2 *( x ) 



(8) 



The final transformation which eliminates the need 
for an unphysical time variable follows similarly to the 
Hamiltonian for the Gaussian thermostat . First add 
a constant to <f> so that the initial (and hence at all times) 
value of Hu is zero. It is easily verified that multiplying 
a zero Hamiltonian by an arbitrary function scales the 
time, but has no other effect on the equations of motion. 
Thus the final form of the Hamiltonian is 

H F (X; P; t) = ^-e*W + 0(X) e -*^ . (9) 

This transformation is equivalent to multiplying Hn in 
Eq. (|l|) by s, which is a particularly simple method of 
generating Eqs. without the use of unphysical time 
variables. 

It is also convenient to introduce a few more 37V + 1 
dimensional vectors, 



V = Pe $ 
F = -V$ 
f = -V<£ 



(10) 

(11) 

(12) 



(13) 



in terms of which the condition Hp = becomes 
V 2 

- + ^° ' 

and Hamilton's equations of motion reduce to 

X = V (14) 

V = V 2 F-F-VV + f . (15) 

As an example, let us consider the Nose-Hoover oscil- 
lator of Ref. Q . For a single particle in a one dimensional 
harmonic oscillator potential, ip — muj 2 q 2 /2. Then the 
Hamiltonian becomes 

H F (q, a; ir,p a ;t) = {it 2 /{2m) + p 2 J{2Q)) e~° 

+ {mu 2 q 2 /2 + gkTa) e a , (16) 

and the vectors we introduced above are 



X = {q^/m, vQhis) 

V = {7T/{ S V^), P a/{sVQ)) = {P/V™~, VQO 

F = (0,1/VQ) 
f = (~LU 2 q\/m,—gkT/i 



Q) 



(17) 

(18) 
(19) 
(20) 



In this form it is straightforward to show that Eqs. 
( |i"^ , p"5|) are equivalent to Hoover's form of the equations, 
Eqs. (|2|j|). The decoupling of the s (or a) equation occurs 
because F and f are independent of s. 
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III. CONJUGATE PAIRING 

The equations of motion given in the previous section 
are now in a form suitable for a proof of the conjugate 
pairing rule for Lyapunov exponents. Note that the coor- 
dinates used are the same (apart from constants related 
to the masses) as the Nose- Hoover thermostat, except 
that s is replaced by a. The components of V are pro- 
portional to the physical momentum p and £. Since the 
equations do not depend on er, the Lyapunov exponents 
obtained contain one extra zero exponent, but otherwise 
are the same as the other 67V + 1 equations considered 
separately. 

The argument given here closely follows Ref. pi]| , so 
the presentation will be correspondingly brief. It is con- 
venient to group X and V together to form a point T in 
67V + 2 dimensional phase space. Time dependent ma- 
trices T and L are defined, giving the infinitesimal and 
finite evolution of linear perturbations ST as follows. 



6t(t) = T(t)ST{t) 
ST(t) = L(t)5T{0) 



(21) 
(22) 



The Lyapunov exponents are defined as the logarithms 
of the eigenvalues of A, where 



A = lim (L T (t)L(t)) 

t—>oc 



l/(2t) 



(23) 



Two of the Lyapunov exponents are zero due to the fact 
that perturbations along the flow simply add a constant 
to the time, and the conservation of the value of the 



Hamiltonian, Eq. (13), which is the total (scaled) energy 
of the system and reservoir. Of course the energy of the 
system alone is not conserved, as it maps out a canonical 
distribution. 

The conjugate pairing clearly does not include these 
two exponents, so perturbations along the flow, and 
which alter the value of the Hamiltonian must be elimi- 
nated before the pairing is apparent. This is achieved by 
considering 67V perturbations, none of which are along 
the flow or alter the total energy. They are measured with 
respect to a basis in a 67V dimensional subspace which 
rotates with the trajectory in order to enforce these prop- 
erties. In particular <5X is always perpendicular to V (in 
37V + 1 dimensional coordinates) , and 5V has a compo- 
nent parallel to V which is fixed by energy conservation. 
As the perturbed trajectory evolves, it will always have 
the same conserved energy, but <5X may not remain per- 
pendicular to V. The perturbed trajectory will, however 
cross the 67V dimensional space at some time t different 
to i, so the above conditions may be enforced by allow- 
ing the perturbed trajectory to evolve at a rate different 
to the original. The Lyapunov exponents obtained using 
this approach are the same as for the full 67V + 2 dimen- 
sional space, with the exception of the two zeros. These 
issues are discussed in more detail in Rcf. ^J, the only 
difference being that here the phase space is not com- 
pact. The arguments follow through exactly the same, 



however, if ergodicity is assumed, as it is when deriving 
the canonical distribution (Sect. |). 

The first step is to choose an orthonormal basis in 37V + 
1 dimensional space, {eo,ei}, with i ranging from 1 to 
37V, as it will henceforth, ep is parallel to V, 



V = Veo 



(24) 



with V = | V| and all the are perpendicular to V. This 
basis is used for both X and V subspaces. The equations 
of motion (jlj-^) determine the time evolution of V and 



V = f • e 

e = Y / (VF + V- 1 f)-e i e i 



(25) 
(26) 



The equations of motion for the other basis vectors are 
somewhat arbitrary, but the natural choice which pre- 
serves the orthonormal character is parallel transport 
along the trajectory, 

e i = -(^F + ^- 1 f)-e i e . (27) 

The perturbed trajectory subject to the above condi- 
tions then becomes 



X' = X + 5X i e i 

i 

V' = V + V' 1 f • ^SXi + 8Vn 



with equations of motion 
dt 

d i fn t i ii i 

— V = V 2 F F VV + f 

dt 



(28) 
(29) 

(30) 
(31) 



Here, F and f are the values at the perturbed positions, 
that is, 



F' = F + SXtViF 

i 

f' =f + ^5JQVif . 



(32) 
(33) 



Substituting Eqs. (H@||,|3|) into Eqs. @|T]), ig- 
noring quadra tic per turbations, simplifying with the help 
of Eqs. (njl^ , B5f|27| ) , and taking components in the di- 
rections of eo and the leads to 67V + 2 equations. One 
of these is not independent of the others, due to energy 
conservation. One relates t and t, and the remaining 67V 
determine the evolution of the perturbations: 



dt 



-^r = 1 + 51 (F + 2V~ 2 f) ■ ei SXi 



(34) 



SXi = SVi (35) 
SVi = SXjej ■ (V 2 VF + Vf - V 2 FF 



3V- 2 ff - Ff - £F) • e, - VF ■ e 6Vi 



(36) 
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Note that VF ■ eo is simply £, and the terms contain- 
ing gradients of forces are symmetric, since F and f are 
derived from potentials. From these equations, the in- 
finitesimal evolution matrix T for the restricted 6-/V di- 
mensional space may be read off as 



T 



I 

M -(I 



(37) 



where each of the elements are 3N x 3iV submatrices. M 
is symmetric because F and f have been derived from a 
potential, and and / are the zero and unit matrices, 
respectively. T satisfies the equation 



T J + JT = -(J , 
where J is given by 



J = 



I 

-I 



(38) 



(39) 



From this point the analysis is exactly the same as 
Ref. |21|. Eq. ( |38| ) leads to similar relations for L and 
L T L, and hence the eigenvalues of A. The end result is 
that for each exponent A, C — A is also an exponent, with 



C 



<(>t 



(40) 



that is, minus the time average of £. Thus the conjugate 
pairing rule holds in the (q, a, p, £) variables, with the 
exception of the two zero exponents. As noted before, 
none of the equations of motion depend on a, so that, 
omitting the a equation gives the same exponents (which 
pair, as shown above), with only one zero exponent. 

Conjugate pairing has now been shown for Gaussian 
and Nose-Hoover thermostats which act on the kinetic 
energy. Numerical simulations [^4| suggest that the con- 
jugate pairing rule holds also for Gaussian thermostats 
which keep the internal energy (rather than the kinetic 
energy) constant. 
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